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Abstract 

For a model of nonlinear elastodynamics, we construct a finite volume scheme which 
is able to capture nonclassical shocks (also called undercompressive shocks). Those 
shocks verify an entropy inequality but are not admissible in the sense of Liu. They 
verify a kinetic relation which describes the jump, and keeps an information on the 
equilibrium between a vanishing dispersion and a vanishing diffusion. The scheme pre¬ 
sented here is by construction exact when the initial data is an isolated nonclassical 
shock. In general, it does not introduce any diffusion near shocks, and hence nonclas¬ 
sical solutions are correctly approximated. The method is fully conservative and does 
not use any shock-tracking mesh. This approach is tested and validated on several test 
cases. In particular, as the nonclassical shocks are not diffused at all, it is possible to 
obtain large time asymptotics. 


Key words and phrases: nonclassical shocks, undercompressive shocks, kinetic relation, 
finite volume schemes 
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Introduction 

Hyperbolic systems of conservation laws often arise as limits of systems of partial differential 
equations including small scales effects, like diffusion and dispersion, when the parameters 
driving the small scales effects vanish. For example, the compressible Euler equations are 
the (formal) limit of the compressible Navier-Stokes equations, when the diffusion param¬ 
eter tends to zero. However, the solutions of the hyperbolic system do not inherit all the 
properties of the augmented system, and it particular they are not always smooth. More 
importantly, in the class of weak solutions, the Cauchy problem for the hyperbolic system 
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may admit infinitely many solutions, plenty of them being non relevant toward the aug¬ 
mented system. Thus it is necessary to add some criterions that select the solutions that 
are indeed limits, as the parameters vanish, of the solutions of the augmented diffusive- 
dispersive system. One of this criterion yields to the selection of classical solutions, which 
roughly speaking corresponds to the case where diffusion is dominant. But when two small 
scale effects are of comparable strength (say diffusion and dispersion) nonclassical solutions 
appear. They contain shocks that do not verify the classical criterion but do arise as limits 
of the diffusive-dispersive system. 

Let us sketch a general framework for nonclassical solutions of hyperbolic systems. Con¬ 
sider the Cauchy problem for a system of conservation laws 

(dtU{t,x) +d,^f{U{t,x)) =0, 

\u{t = 0,x) = U%x), 

where t is positive, x belongs to M, U is the vector regrouping the > 1 unknowns belonging 
to a open convex set Q of and / = > R^ is the flux function, which we assume 

to be smooth. We suppose that System ([T]) is strictly hyperbolic, i.e. that the Jacobian 
matrix of / is diagonalizable with simple eigenvalues Xi{U) < ■ ■ ■ < Xn{U). We denote by 
ri(U) a right eigenvector associated to the eigenvalue \i{U). 

Solutions of hyperbolic systems are in general not smooth, even when the initial data is 
smooth. Discontinuities appear in hnite time and it is necessary to consider weak solutions. 
Doing so, uniqueness is lost, and an additional criterion has to be imposed to select a unique 
solution. For example, the vector of unknown U is asked to verify, in addition to the initial 
System m, a so called entropy inequality 

dtU{u) + d,w{u) < 0. (2) 

The entropy U and the entropy flux W usually come from considerations on an augmented 
version of the system under study, where small physical terms like diffusion and dispersion 
are not neglected. Such an augmented system has additional properties, for example the 
total energy is conserved. Inequality ([2]) states that at the limit, this energy is dissipated. 
In general, from the augmented system only one entropy inequality is deduced, thus the 
solutions of © are not expected to verify other entropy inequalities than ([2]). When the 
system is genuinely nonlinear, i.e. when 

Vf ,n}, VC/gD, S/Xi{U)-niU)^0, (3) 

and when the entropy U is stricly convex, the entropy inequality ([2]) contains enough infor¬ 
mation to select only one weak solution of ([T]). In the sequel we suppose that the entropy 
U is strictly convex. In the general case where ([3]) does not hold, the addition of a single 
entropy inequality is not sufficient to select a unique weak solution. 

The question of uniqueness of the weak solution is strongly related to the choice of 
a criterion to decide wether a shock is admissible or not. Two states Ul and Ur of D 
are linked by a shock if there exists a real s, called the speed of the shock, so that the 
Rankine-Hugoniot relations 

s{Ul - Ur) = f{UL) - f{UR) (4) 
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hold. The shock is entropy satisfying if it verifies ([2]) in a weak sense, i.e. if 

s{U{Ul) - U{Ur)) - {W{Ul) - W{Ur)) < 0. (5) 

When the system is genuinely nonlinear (i.e. when (|3)) holds), ([5]) is equivalent to the 
Lax ineoualities [Lax m 

e {!,■ ■ ■ ,n} ■. \i{UR) < s < \i{UL)-, s < Xi+i{UR) and s>Xi-i{UL), (6) 

which themselves are an extension to systems [N > 1) of Oleinik’s criterion |01e57| . 

In |Liu74| . Liu generalized this criterion to the non genuinely nonlinear case where 
VAj • rj may vanish. Liu proved that his criterion selects a unique selfsimilar solution to 
the Riemann problem ([TH2]). We recall that a Riemann problem is the case of a piecewise 
constant initial data 


U\x) = Ul1.<o + UrI^^o, {Ul, Ur) G 

For systems that are not genuinely nonlinear, the Liu criterion is stronger than ([5]) (whereas 
it is equivalent for genuinely nonlinear systems). An admissible shock in the sense Liu 
always verifies ([5]), but some shocks may verify ([5|) without satisfying the Liu criterion. 
These shocks are called nonclassical. Moreover, unlike in the genuinely nonlinear case, the 
addition of the single entropy inequality ([2]) is not enough to regain uniqueness of the weak 
solution, because too many shocks are allowed. In other words, the inequality entropy ([2]) 
does not contain enough information on the small scale effects in the augmented system. 
Uniqueness can be regain by strongly constraining the nonclassical shocks to verify an 
algebraic relation of the form 

Ul = ^^{Ur). 

The function is called a kinetic relation. 

It has been proved that the addition of such a kinetic relation selects a unique weak 
solution to the Riemann problem for several models, including models of phase transitions 
in solids in |AK91j . |LTni| and |LTn2] or in liquids [LTDO] . [HLOOj, |SY95] and a model of 
magnetohydrodynamics [HLOO] . The link between kinetic relations and augmented diffusive 
dispersive equations is explored in |,IMS95] . |LeFn2| and |HL97) in the scalar case and 
in |Sle89] for gas dynamics with a Van der Waals pressure law. The literature on the 
subject is so wide that we only cited a few key references. The interested reader will find a 
more comprehensive bibliography on this topic in |LeF02| . 

The numerical approximation of nonclassical solutions is challenging, because they are 
by nature very sensitive to the equilibrium between the small scale effects. The stake is 
to preserve this equilibrium at the numerical level. Roughly speaking, usual finite vol¬ 
ume schemes introduce a numerical viscosity which destroys the equilibrium between the 
small scales effects. The diffusion becomes dominant and the schemes converge toward the 
classical solution, even though they are based on nonclassical Riemann solvers. 

To overcome this difficulty, two types of schemes have been proposed. On the one 
hand, it is possible to use high order schemes consistent with the augmented system (see 
for example |HL98) . |LMR02] . |LR00| and |KR10| L As outlined in |HL98| , it is necessary 
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to discretize the flux of the hyperbolic part with a high enough order, once again to keep 
the exact balance between the small scales. On the other hand, some schemes relies on 
the tracking of nonclassical shocks (see for example [CL03) . |CG08| . |BCLL08] . |Perll| . 
|OORR12| and |0DM014) ). In that case the kinetic relation is taken into account the 
scheme, typically through the use of an exact nonclassical Riemann solver. 

The conservative finite volume scheme presented here belongs to this last category. It 
is built to be exact when the initial data is a nonclassical shock. In general, it does not 
introduce any numerical diffusion near nonclassical shocks (a particular class of nonclassical 
solutions) and turns out to capture correctly nonclassical solutions. This scheme extends to 
hyperbolic systems the discontinuous reconstruction scheme introduced in the scalar case 
in |BCLL08) and recently used in |CDMG14] . 

In the first section, we present the hyperbolic system admitting nonclassical solution for 
which we construct the scheme, namely the model of nonlinear elasticity studied in [LT01| 
and |HL00| . The second section is devoted to the construction of the scheme itself. We 
present a way to select cells in which a nonclassical shock is reconstructed, explain how 
this shock is reconstructed and give the numerical flux. Eventually, we prove that the 
scheme is exact when the initial data is a nonclassical shock. In the third and last section, 
we propose several test cases involving nonclassical shocks. It appears that the proposed 
scheme capture nonclassical shocks very sharply. Let us outline that the scheme is fully 
conservative (unlike Glimm type schemes |CL03) 1 and uses a fixed grid, which allows the 
computation of solutions containing several interacting nonclassical shocks. 

Acknowledgment The author warmly thanks Frederic Lagoutiere for his support and 
advice. 


1 The Riemann problem for a nonlinear elasticity model 


This paper is devoted to the numerical approximation of the solutions of the system of 
conservation laws 

✓ 

dtv - dx(T{w) = 0 , 
dtw — dxV = 0, 
v{t = 0, x) = v^{x), 
w{t = 0,x) = lo^(x), 


(7) 


where the stress cr is twice differentiable, and verifies 


wa"(w) > 0, (t'{w) > 0 and lim a' {w) = +oo. 

|io|—>-+oo 


( 8 ) 


We are interested in weak solutions of ([7]) that are morally limits when e tends to 0^ of the 
augmented system 

j dtv'^ - dx(T{w^) = edxxV + ae^dxxxW, 

I dtw - dxV = 0. 


(9) 
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Figure 1: Notations used for solving the Riemann problem. is the first Liu shock 

reachable from w. The entropy dissipation between w and is zero. The states in red 
can be linked to rc by a Liu shock, the blue ones by a nonclassical shock. 


The parameter a is positive. From Q we deduce the following entropy inequality for 
System ([7]) (see |LeF02| L 

dt + J a{z)dz^ + dx {-va{w)) < 0. (10) 

In this section we briefly recall the results of Thanh and LeFloch [LTOl) on the Riemann 
problem for System ([7]), i.e. the case where 


r;0(a;) = ULla,<o + Vr1x>0, 

w^{x) = Wl1x<Q + Wr1x>Q. 


( 11 ) 


We adopt the notation of this paper. It is recalled on Figured] 

If the left state {vl,wl) and the right state {vr,wr) are linked by a 1-shock, its speed 
is equal to —s{wl,wr), where 


s{wl,wr) 


^ ajwL) - crjwR) 

WL-WR 


( 12 ) 


and vr is given by 


VR = VL + s{wl,wr){wr - wl) ■■= T-Li{wr,vl,wl). (13) 
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If they are linked by a 2-shocks, its speed is +s{wL,wii), and 

VL = VR- s{wl,wr){wl - wr) := 'H2{wl,Vr,Wr). 

Graphically, the speed of a shock s{wl, wr) corresponds to the slope of the segment joining 
to {wR,a{wR)). 

Let us now focus on when a 1-shock is admissible or not, either in the sense of the 
entropy dissipation m or in the sense of Liu. The results for 2-shocks are the same, with 
the roles of wl and wr reversed in what follows. A careful study of Equation (IlOp applied 
to the shock yields that it verifies the entropy inequality if and only if 

wl < WRWR or WRWL < 

where is the real having the opposite sign than wr for which the entropy dissipa¬ 

tion jb]) vanishes. The function is its own inverse. 

For System ([7]), a shock verihes the Liu criterion if and only if for all w between wr and 

WR, 

a{w) - a{wR) ^ a{wR) - a{wR) 

W -WR ~ WR-WR 

Geometrically, it means that the segment joining {wR,a{wR)) to {wR,a{wR)) is above the 
graph of a if wr > wr, and below if wr < wr, see Figure [TJ It follows that there exists an 
invertible function <I>^ such that a shock verihes the Liu criterion if and only if 

w\ < wrWr or wrWr < ^~'^{wr)wr. 

The fact that System ([7]) is not genuinely nonlinear yields 

wr^~'^{wr) < wr^'^^{wr). 

Thus, if Wr lies in between ^1^{wr) and ^~\wr), the shock between {vr,wr) and {vr, wr) 
dissipates the entropy m but is not admissible in the sense of Liu. Such shocks are called 
nonclassical shocks. 

Many Riemann problems for System dZHinD admit an inhnity of solutions (see [LTm] . 
Section 3). Uniqueness can be obtained by imposing a kinetic relation that strongly con¬ 
strain the nonclassical shocks by imposing 

{ Wr = ^^’^{wr) for a 1-nonclassical shock, 

Wr = ^^’’^(wr) for a 2-nonclassical shock, 

where et both verify 

< w^^’'^{w) < w^\w) 

but are not necessarily equal. We now state the main result of |LT01| . 

Theorem 1.1. [Thanh, LeFloch] If the kinetic functions and are monotone de¬ 
creasing, the Riemann Problem (fTtliOlfiTlj has a unique selfsimilar solution. 
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The exact Riemann solver associated to Theorem 11.11 is the foundation of the scheme 
presented in the next section. Moreover, we will use the inner structure of the Riemann 
solution. Thus for the sake of completeness, we describe below the forward 1-wave and the 
backward 2-wave. 

Proposition 1.2. [Thanh, LeFlochj Let {vl,wl) be a fixed left state, and [vm^wm) be a 
state linked to {vl,wl) by a 1-wave. Then the nature of this 1-wave is the following: 

• ifwLWM > it is a classical shock; 

• if 0 < wlWm < w\, it is a rarefaction wave; 

• if < wlWm < 0; it first contains a 1-rarefaction linking {vl,wl) to 

the state which is itself linked to {vm,wm) by a 

1- nonclassical shock; 

• if WlWm < 'WlT~^’^{wl), two cases arise: 

— if —s{wL,(p'^’^iwM)) < —s{^P^’^{wm)),'Wm), it first contains a 1-classical shock 
linking {vl,wl) to the state {'Hi{i.p^'^{w),vm-,wm)-,T^’^{wm)), which is itself 
linked to {vm,wm) by a 1-nonclassical shock; 

— otherwise, it is just a 1-c/asszca/ shock (in which w changes sign). 

Proposition 1.3. [Thanh, LeFloch] Let {vRjWLi) be a fixed right state, and {vm,wm) be 
a state such that (vm^wm) CLnd {vrjWr) are linked by a 2-wave. Then the nature of this 
2-wave is the following: 

• ifwRWM > it is a classical shock; 

• if 0 < wrWm < it is a rarefaction wave; 

• ifwii(p~^’'^{wii) < wrWm < 0, it first contains a 2-nonclassical shock linking {vm, wm) 

to the state which is itself linked to {vii,wii) by a 

2- rarefaction wave; 

• if wrWm < WLii.p~''’'‘^{wR), two cases arise: 

— if s{wmt‘^'‘^{wm)) < s{ip^’'^{wM)),WLt), it first contains a 2-nonclassical shock 
linking {vm,wm) to the state {(H 2 {t^’‘^{w),vmiWm)^T^'^{wm)), which is itself 
linked to (vr, wr) by a 2-classical shock; 

— otherwise, it is just a 2-classical shock (in which w changes sign). 

The distinction between the last two cases in the last point of the enumerations of 
Proposition 11.21 and 11.31 can be reexpressed geometrically, as illustrated on Figure [2l 

Proposition 1.4. The exists a function (p^ such that for i G {1,2}, for all real w and 
for all real wm such that wwm < w(p~'^’'^{w), the quantity \s{w,(p'^’'^{wm))\ is larger than 
\s{(p'^’^(wm)),wm)\ if ond only if wip^’'^(wM) is larger than wip'^(w,WM)- In, particular if 
0 > wwm > wmt'^{w,wm), the solution is a classical shock (in which w changes sign). 
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Figure 2: The solution contains a 2-nonclassical shock (in blue) followed by 2-classical shock 
(in red) if their speeds are correctly ordered, i.e if the slope of the blue segment is smaller 
than the red one (right). This depends on the relative position of and ,wl)- 


2 The scheme 


We now describe the scheme. The time discretization is denoted by = 0 < < • • • < 

t” < • • •, and at each time t"", the space is discretized with cells having all the same size 
Ax. We denote by their extremities and by x” there centers. The centers of the cells 

will move at each time step but every cell will always be of size Ax. The speed of the mesh 

is defined by ■ Integrating the first line of ([7]) over the set 


{{t, x),e < t < + Kesht <X< + ''^meshO 

we obtain the integral formula 


i-i/2 


f^j + l/2 , f^j + 1/2 

/ v{t'^~^^,x)dx = / v{t‘^,x)dx 

Jx" , 

- / a:^+l/2 + 'KSesh(^ 

- - lCe,h^)(s, x]_^/2 + ''Cesh(s ’ 


- r))ds 
r))ds 


The formula is identical for the variable w with the flux {—v — I^esh^) instead of {—a{w) — 
I^esh^- principle of finite volume methods is based on this integral formula. A finite 





volume scheme writes 




_ „ „ / rn,v rn,v \ 

~'^j~ Ax yj j+l/ 2 ~ J j-1/2)^ 
(_pn,w _pn,w \ 

Ax w j+1/2 “ Jj-112)^ 


= % - 


where plays the role of the mean value 


1 r7+i/2 

— / vir,x)d3 
Ax ^ ^ 

^j- 1/2 


and f'^^i /2 pl^-ys the role of the flux on the edge x = Xj_ij 2 + l^esh(® “ 
1 


fn+l _ 


[ i-cr{w) - V^esh^)('S, a^”+i/2 + K^esh('S “ t''))ds 


(15) 


(and similarly for tc” and f^^^ 2 )- scheme is initialized with the exact average of the 
initial data 

n 1 n n 1 /'^i+1/2 n 

Vj G Z, / v^{x)dx and vJ^ = / w^{x)dx. (16) 

A0_i/2 Jx0_j/2 


A particular finite volume scheme is characterized by a particular choice of a formula ex¬ 
pressing the numerical fluxes f^^i /2 f ^^/2 ™ terms of the (u(()fcez and {w'^)k£Z- Our 

choice is described in the next sections. 

Recall that our aim is to derive a scheme which is exact when the initial data is a 
nonclassical shock. In that case, the initial sampling (1161) introduces a small amount of 
dissipation (unless if by miracle the shock initially falls on an interface). The numerical 
initial data contains an intermediate value that does not correspond to any pointwise value of 
(u*^, tc^). More generally, at the end of every time step, finite volume schemes contains a L^- 
projection on the mesh. Details of the solution are lost in that step, and intermediate values 
are created in the shocks profile. The main idea of the discontinuous reconstruction scheme 
is to rebuilt, from the mean values (u”,rc”)jgZ) an initial data that contains more details. 
This idea is also the foundation of well-known other schemes, like the MUSCL scheme [vL97| 
and its central version |NT90| . Our reconstruction consists in adding nonclassical shocks in 
some cells in which we are able to detect that a nonclassical shock where lying inside the 
cell before the L^-projection. Thus the reconstruction will be very precise near nonclassical 
shocks, and not on the smooth parts of the solution, taking the opposite strategy of [vL97| 
and |NT90| 

Following this idea, our scheme can be decomposed in three elementary steps. 


• First, we detect the special cells which are associated with nonclassical shocks. This 
detection step is described in Section [2.11 


• Then, we reconstruct nonclassical shocks inside those particular cells, in the sense 
that we replace the mean values {v^,w'j) by piecewise constant functions having the 
form of a nonclassical shock. This procedure is explained in Section [2.21 
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• Eventually, the numerical fluxes are computed by letting the reconstructed nonclas- 
sical shocks evolve during the time step — t”. The use of a moving mesh makes 
this computation easy, see Section [2.31 

In the last section, we extend the scheme to reconstruct classical shocks as well. 


2.1 Detection of nonclassical shocks 


The key idea of the scheme is to see, whenever it is possible, each mean value {v^,w'j) 
produced by the finite volumes scheme as the average of some nonclassical shock located 
somewhere inside the cell. Of course, where and are linked by 

a nonclassical shock, it is the one that should be reconstructed. In general, 
and are not linked by a single nonclassical shock, but by the full pattern of 

two waves, each of them likely to contain a succession of a classical shock or a rarefaction 
followed by a nonclassical shock (see Propositions 11.21 and II.3p . In that case, we reconstruct 
one of the nonclassical shock appearing in the solution of the Riemann problem. It is chosen 
thanks to the following Lemma. 

Proposition 2.1. If the states {vl^wl) and {vR,wji) are linked by a 1-nonclassical shock, 
then 

WLWR < 0 , {wL - wr){vl -vr) >0 and wrwr > WRip^{wL,WR). (17) 

If they are linked by a 2-nonclassical shock, then 

WLWR < 0 {wL - wr){vl - vr) <D and wrwr > wl‘I>^{wl-,wr). 


Proof. The condition wlWr < 0 always holds through a nonclassical shock. The second 
one is a straightforward consequence of the Rankine-Hugoniot relations (H]). In the case of 
nonlinear elastodynamics ([7]), they write 


{-iys{wL,WR){vL - vr) = a{wR) - a{wL)] 
{-iys{wL,WR){wL - Wr) =Vr- vr. 


(18) 


with i = 1 for 1-shocks and i = 2 for 2-shocks. The positive real s{wl,wr) is defined in 
Formula (fT2ll . Eventually, the last condition is a reminder of Proposition 11.41 □ 


Definition 2.2. For all integer j and all positive integer n, we denote by {vj'^,wy^) the 
intermediate state appearing in the Riemann problem with left state rc4_^) and right 

state {Vj‘_^i,Wj'_^i). The left and right desired reconstructed states, denoted respectively by 
{v^l^'^^^l) defined as follows. 


If wy_^ 


^i+i 

n 


< 


0, {wy_^ - wy^^)ivy_^ - vy^^) > o. 


Wj-iwy+i > 


w 


j+i 


(p'^{wy_ 




and wy_iwy_,^ < 0, the solution contains a 1-shock in which w changes sign. 


— If this shock is nonclassical, it is preceded by a classical wave and we set 






v^ w"‘ '' 

’ J,*’ 




{v 




) =i< 


4 yj"- ) 
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— If this shock is classical, it is the only contribution in the 1-wave and we set 


(■ 








If < 0, (w^_, - < 0, 

and < 0, the solution contains a 2-shock in which w changes sign. 

— If this shock is nonclassical, it is followed by a classical wave and we set 




— If this shock is classical, it is the only contribution in the 2-wave and we set 

=K+i,<+i)- 

In the other cases, we do not detect any relevant nonclassical shock and we set 


Kr^^Ir) 


) =( 


= (Vj,Wj 
= 


Remark 2.3. If there is a nonclassical shock in the Riemann problem with left state 
and right state then the first two three conditions holds by Proposition 12.11 

Note that they are numerically very easy to check. The last condition insures that the 
Riemann problem indeed contains the expected shock; we verify it only on the cells that 
passed the hrst three tests. 


2.2 Reconstruction 

Once a nonclassical shock has been detected, it is placed inside its cell by conservation of 
the variables v and w. We reverse the averaging step of finite volume schemes by replacing 
the mean value u” by the piecewise constant function 

'^rec(^) ~ ^j,L'^x<Xj_i/2+d'J’'" ^j,R^x>Xj_i/2+d”’'' 


with 




= Ax 




(19) 


If belongs to (0, Ax) we have 


f 


Kecix) dx 


Axvj 
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and no mass is loss when replacing the mean value by u”c inside the j-th cell. 
Reasoning similarly for the w variable, we replace by 


Kecix) = wIl^ 




with 


= Ax 


'^1,R 


(20) 


We now state one trivial but crucial property of this reconstruction procedure. 


Proposition 2.4. If and are linked by a single nonclassical shock, 

and if it exists a in (0,1) such that 


(u7,0 = + (1 - a)(u;+i,<+i). 


then 


{vIl, wIl) = i^lR^ <r) = <+i) d^ = dj’" = aAx. 

In particular if the initial data is a nonclassical shock, we have uj?g(. = u® and = 
w^: the numerical diffusion introduced by the initial sampling (|16p is cancelled by the 
reconstruction. 

In general the two distances and d^’^ are different, and it is possible that at least 
one of these distances does not belong to (0, Ax). In that case we cancel the reconstruction, 
considering that seeing (u”, tc”) as the mean value of the detected nonclassical shock is not 
relevant. 

Definition 2.5. The left and right reconstructed states are defined by: 


yjn if dY € (0, Ax) and € (0, Ax), 

otherwise. 


( 21 ) 


( n n N I if dY e (0, Ax) and dj’"' G (0, Ax), 

^ \{vf,wY) otherwise. 


( 22 ) 


According to the Rankine-Hugoniot relations ([3]), the nonclassical shock reconstructed 
in cell j has velocity 


„n _ ^j,R '^j,L 
Sj - „ 




and we set arbitrary s^ to a'{wj) when no reconstruction is performed (which reads 


Yl = Yr = Y") 
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2.3 Advection of the reconstructed discontinuities 


The fluxes are computed by letting the reconstructed nonclassical shocks evolve during 
the time step and by computing exactly what goes through the interfaces. However, two 
discontinuities reconstructed in adjacent cells can interact and the waves resulting from the 
interaction can meet the line x = within the time step. It follows that if we want 

to use a fixed grid = x'j, or equivalently = 0), the flux along the interface 

X = cannot be computed without resolving the wave interaction, which is obviously 

extremely costly. 

This can been avoiding by using a moving mesh. Let us recall that by moving mesh, we 
mean that the centers of the cells are moving from time to time, but that their size remains 
constant equals to Ax. Thus the numerical difficulties of handling cells with different and 
varying widths is avoided. The mesh speed is chosen such that it is larger than the maximum 
of the waves speed: 

iKSeshI > Kaves = max (23) 


and the time step such that a wave cannot cross more than an entire cell during the time 
step: 


^n+l _ < 


Ax 


IK 


mesh 


' ' in 


(24) 


These two hypothesis insure that any wave created at time t"' inside the j-th cell can 
only cross the right interface x = “1“ Knesh^ Knesh < ^ interface 

X = + Knesh^ Knesh > 0. It also imply that if two waves interact inside the j- 

cell during the time step, the waves resulting from their interaction will not have time to 

catch up the left or the right interfaces. This is depicted on Figure [3j Thus, computing 
the flux along the j + 1/2-th interface x = + Knesh^ boils down to compute the time 

at which the nonclassical shock reconstructed in cell j crosses the interface if Kiesh ^ 
(in cell j + 1 otherwise). On Figure [3l we also see that if no reconstruction is performed 

in cells j and j + 1, the flux is trivial to compute and is given by a simple left or right 

decentering (depending on the sign of Knesh)’ which exactly corresponds to the staggered 
Lax-Friedrichs flux. 

We now compute the flux in details, when l^iesh negative and verifies (|23p . like in 
Figure [3j Then the nonclassical shock reconstructed in cell j may only cross the right 
interface x = a^Ki/2 Knesh(^ ~ discontinuity is not located at the same place for 

the variables v and w, therefore we compute two different crossing times 


T 


Ax — d 


\n^w 

i-i 


vj^n _ 

j+l/2 ~ gU yn 

J — 1 mesh 


and 


rpW,n 

i+1/2 


Ax - 


A-i 


-K 


if K 


mesh 


< 0. 


mesh 


Once again. Condition (I23p insures that the flux passing through the j + l/2 interface only 
comes from waves arriving from the cell j. It can be clearly seen on Figure [3] that it is 
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^n +1 



Figure 3: Computing the flux through the dashed interfaces can be done without solving 
waves interactions (interface j + 1 / 2 ), and is trivial when no reconstruction is performed 
(cell j - 1 ) 


piecewise constant, given by 


' rn,V _ 

•/j+ 1/2 - 


nn^W _ 


- VZs^vlLm- - min(AtM;;’i/2)) 
(-vIr - min(Ar, 

H-vIl - " min(Ar, 


if ^mesh < 0 (25) 


where we denote by At”’ the n-th time step — t”. 

When V/(gsjj is negative, the reasoning is the same but the shock reconstructed in the 
j-th cell now crosses the left interface x = f/nesh(f “ crossing times are 

given by 


rpW,n 

^i- 1/2 


_ „n 
mesh '^j 


and 


rpW,n 

^i- 1/2 


K 


mesh 


— S. 


if y” 


mesh 


> 0 . 


and the fluxes are 


pn,v _ 

Jj-ll2 


nn,W 

Jj-1/2 


(-^Kl) - h2esh^/:L)nfln(At”,7;_:’;/2) 
i-y^L - Kes^y^lL) min(At”, i;_:7/,) 

H-vIr - - min(At”,i;_:-2)) 


if > 0 (26) 


Remark 2.6. When no reconstruction is performed, the fluxes (|26p and (|25H coincide with 
the staggered Lax-Friedrichs fluxes, i.e., to a simple left or right decentering depending on 
the sign of The values assigned to the distances and and to the speed s” 

(and thus to the crossing times ^^ 4 .^/ 2 ) °f importance. 
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Remark 2.7. The idea of using a moving mesh to simplify the computation of the fluxes 
goes back to the Lax-Friedrichs scheme. The Rusanov scheme follows the same idea by 
using a local mesh speed to make the scheme less diffusive. Higher order extensions 

based on piecewise polynomial reconstructions are possible, see for example |NT90| . 

2.4 Detection of classical shocks 

It is also possible to detect and reconstruct shocks in which w does not change sign. Those 
shocks are always classical. Their detection is based on the following proposition. 

Proposition 2.8. If the states {vl,wl) and {vrjWr) are linked by a 1-shock in which w 
does not change sign, then either 

wr < wr < 0 and vr < vr 


or 


0 < Wr < Wr and Vr > Vr. 

If they are linked by a 2-shock in which w does not change sign, then either 


Wr < Wr < 0 and vr < vr 


or 

0 < Wr < Wr and Vr > Vr. 

There is no conflict with the detection of nonclassical shocks of Proposition 12.1[ and we 
can straightforwardly extend Definitions 12.21 and 12.51 to take into account those shocks. 


3 Exact approximation of isolated nonclassical shocks 

The aim of this section is to prove that the scheme described above is exact when the 
initial data is an isolated nonclassical shock, i.e. m with left state {vr,wr) and right 
state {vr, Wr) verifying the Rankine-Hugoniot relation (1181) and constrained by the kinetic 
relation (|14p . We recall that in that case the exact solution is 

{ V (t,x) '^Llx<{—iys{wL,wii)t '1~ s(wi:^,wii)t 

W (t,x) WRl-x>{—iys{wi:^,wii)t 

with i = 1 for a 1-shock and i = 2 for a 2-shock. 

Theorem 3.1. The scheme described in the previous section is exact when the initial data 
is a single nonclassical shock. In other words, the numerical solution is the -projection 
on the mesh of the exact solution at time t”.’ for all n > 0 and for all j € Z, 

1 rl+i/2 1 /'^"+i/2 

< = — / v^^^{e,x)dx and wf = — w^^^ir,x)dx. 

I I 
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w = 

0 

Wq 


w = 0 



WR 

-1 

0 

1 

WR 




Figure 4: Structure of w when a 1-phase transition is detected in cell —1. In that case, the 
Riemann solution (on the right) contains three shocks. 


Proof. We prove the result for a 1-nonclassical shock such that wl < 0 < wji, the other 
cases being exactly similar. With the initial sampling m, the property holds true at 
time t^. Suppose that is holds true at some time t”. At this time, we renumbered the 
cells in a way that the discontinuity lies inside the cell numbered 0, and we denote by <5 its 


v-=< 


- 1 / 2 ' have 











VR 

if 

j 

A 

o 



WR 


if 4 

< 

0, 

Lvr + ^iJvR 

if 

j 

= 0, 

and 

w^=< 


Ax—S^,, 

Ax 

if j 

= 

0, 

[vr 

if 

j 

> 1, 



WR 


if j 

> 

1. 

Wr < Wr and Vr 

> 

VR 

, we 

have 

W^l < 

Wq < Wf 

and 

< < 

< 



The 

detection step of Proposition 12.II detects a 1-nonclassical shock in cell 0. By Proposition 12.41 
the nonclassical shock is reconstructed in that cell, and it is be placed exactly at the right 
position in the reconstruction step: 


= {vL,Wl), {Vo^r, wIr) = {VR, wr) 


and 




If no reconstruction is performed in the other cells, our scheme gives the correct solution 
at time by construction. This is easy to check, but a little tedious, and we refer the 
reader to |Agul4| for a detailed proof. 

Suppose that Wq > 0. Then by Proposition 12.81 no classical shock is detected in cell 
1, and by Proposition 12.11 a 1-nonclassical shock is detected in cell —1. Let us focus on 
the solution Riemann problem between the state {vl,wl) in cell —2 and the state (vq,Wq) 
in cell 0. The notation are recalled on Figure HI If w^i > wl, both and 

are larger than wr. Thus it is be impossible to have € (0, Ax) and no reconstruction 
occurs in cell —1. Suppose now that wr \s larger than Then wr and are 

linked by a 1-classical shock, and is linked to ^ by a 1-nonclassical shock. It is 

is impossible to have in (0,1). Indeed, the kinetic function decreases, thus 


wfi^R = (j) > 0 ^’^(wr) = WR 

and there is a classical shock between jj, wfi jf) and {vr, wr). The states wf^ jf) 

and {vr,wr) are both on the 1-wave curve Wf (vr,wr). Theorem 4.5 of |LT01] states that 
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along this curve, v is an increasing function of w. Therefore ^ is larger than vn- On 
the other hand, the Rankine-Hugoniot relations applied to the 2-shock yields that is 

smaller than which contradicts the fact that < vr. 

The case where a 1-nonclassical shock is detected in cell 1 is simpler. We have Wq < 0 
and nothing is is detected on cell —1. Moreover, thus it is not possible to 

reconstruct tc in a conservative manner if w^r < wr. On the other hand, if Wir > wr, the 
solution contains a 2-classical shock and dH yields that v '^r < vr. The Rankine-Hugoniot 
applied to the 1-nonclassical shock implies that < v'^ r- Thus it is impossible that 
belongs to (0, Ax). □ 

Remark 3.2. It is crucial to require that both v and w are reconstructed in a conservative 
manner. In the case where the nonclassical shock is detected on cell 1, it is clear that the 
proof collapses if we authorize outside (0, Ax). Moreover if we authorized d^'^ outside 
(0, Ax), the solution of the Riemann problem between (u”, re”) and {vr, wr) might contain 
a 1-classical shock, a 1-nonclassical shock and 2-rarefaction wave, in which case 

vIl <vr< vIr 

and a reconstruction is performed in cell 1. 

4 Numerical Simulations 

For all the numerical simulations presented below, the stress function is a{w) = + mw 

with m > 0. In that case, 

=—^w, d>^(ui) = —tc and =—w — w'. 

For the kinetic functions we take = d>^’^(t(;) = —/dw, where /3 belongs to [—1/2,1]. 

The case /3 = 1/2 corresponds to the classical solutions ; the choice /3 = 1 corresponds to 
the case where the entropy dissipation (IlOp is zero across nonclassical choice. It does not 
fall in the theory of |LeFr)2| . but the Riemann problem can be solved (see |LTni) i and it is 
possible to explore that case numerically. 

Test 1: Isolated nonclassical shock 

This test case illustrates Theorem 13.11 The initial data is the Riemann problem: 

I v0(x) = -10 * la.<o -b 110 * la;>0, 

|r(;°(x) = -6 * l2,<o -b 9 * lcc>o, 

With m = 1 and /3 = 2/3, the solution is an isolated nonclassical 1-shock. On the top of 
Figure [5l we plot the exact nonclassical solution and the solution given by the reconstruction 
scheme. As expected they are exactly the same. On the bottom of Figure El we plot the 
solution given by the Godunov scheme based on an exact nonclassical Riemann solver. It 
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does not capture the nonclassical solution but the classical one, which in that case in a 
rarefaction followed by a shock. This is a general phenomenon: usual finite volume schemes 
are not able to capture nonclassical solutions. Thus in the sequel we used the Glimm scheme 
to compare our scheme with. The CFL number is set to 0.45, the final time is T = 0.038 
and the space interval [—0.5, 0.5] is discretized with 200 cells. 




10 

8 

6 

4 
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0 
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-4 
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- 8 ^ 


-0.45 -0.4 -0.35 

X 


120r 
100 ■ 


— nonclassical exact solution 
—^reconstruction scheme 


40 - 
20 - 


-0.35 

X 



-0.3 -0.25 



Figure 5: The nonclassical Godunov scheme is unable to capture nonclassical shocks. On 
the contrary, the reconstruction scheme captures isolated nonclassical shocks exactly. 


Interlude: The Glimm scheme 

In the sequel we use the random sampling method of Glimm |Gli65| to compute reference 
solution. The Glimm scheme is built as follow. Let be a sequence of i.i.d. random 

variables uniformly distributed over [0, Ax]. Suppose that at time a piecewise constant 
approximation of the solution is given, and denote by (f, x) e->• x) the exact solution 

at time t > 0 for that initial data. The numerical solution of Glimm at time 
is given by 

Vj e Z, t/”+‘ = (/".(AC,+ r"). 
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It At”’ is small enough, i.e. if it is smaller than the maximum of the waves speed appearing 
in the Riemann problems multiplied by the cell size Ax, is just the juxtaposition of 
Riemann solutions at each interface. 

The main feature of the Glimm scheme is that is does not introduce any numerical 
diffusion, in the sense that the shock profiles are not smeared out at all. This is why this 
scheme has been used in |CL03| to approximate nonclassical solutions. The two drawbacks 
are that the scheme is not conservative and that the exact computation of the solution is 
costly. 

Test 2: A Riemann problem with two nonclassical shocks 

The initial data is now 

(v^(x) = 6 * - 10 * 

|r(;°(x) = 1 * l2,<o + 2 * la;>o. 

The exact solution consists in a 1-classical shock, a 1-nonclassical shock, a 2-nonclassical 
shock and a 2-rarefaction wave. On Figure [6l we plot the exact solution and the numerical 
solution given by the reconstruction schemes. We compare the reconstruction where we 
detect nonclassical shocks only (cf Proposition 12.Ijl . referred to in the legends as RecNC, and 
the reconstruction scheme were the classical shocks are also detected (cf Proposition 12.Sp . 
referred to in the legends as RecNC+C. Both of them capture very sharply the nonclassical 
shocks; the 1-classical shock is much more diffused when it is not reconstructed, and both 
scheme behaves in the same way in the rarefaction wave. The CFL number is 0.45, the 
space interval [—1,1] contains 200 cells. 

Test 3: Perturbation of a classical shock 

This test case is taken from |CL03| . We still have /? = 2/3 but now m = 2. The initial data 
is 

f U°(x) = 1 * lcc <0 -11* la;> 0 , 

|r(;°(x) = (1 e) * lj,<o - 3 * lx>o- 

with e G {0,0.05,0.1} When e = 0, the solution is a 1-classical shock. When e > 0, the 
classical shock is split into a classical shock followed by a nonclassical shock. Their speeds 
are different but close to each other. We plot the solutions at time T = 0.4 on Figure [71 
using 600 cells per unit interval and a CFL number of 0.45. The solutions are very well 
approached when e > 0. When e = 0, a spike appears. This spike corresponds to the state 
linked to {vR,wji) by a nonclassical shock (on Figure [7] we see that it has the exact same 
height as the nonclassical shock appearing for e > 0). Remark that e = 0 is exactly the limit 
between the two cases in the last point of the enumeration of Proposition 11.21 for which the 
speed of the nonclassical and the classical shock coincide (see also Figure [2]). 

Numerically the mechanism is the following. After one iteration in time, an intermediate 
value {vi,Wi) = (avi -t- (1 — a)vji,awL -|- (1 — a)wR) is created. At the second iteration 
in time, the classical shock (in which w changes sign), is perfectly reconstructed in the 
corresponding cell. However, a 1-nonclassical shock is also detected by Proposition 12.11 in 
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Figure 6: A Riemann problem with two nonclassical shocks at time 0.15. 


the cell just before, and this time the reconstruction succeeds. Thus the scheme is not exact 
in that case. Note that the proof of Theorem 13.11 uses the kinetic relation and thus collapses 
when the reconstructed shock does not verify the kinetic relation (and is hence classical). 

This phenomena does not prevent the scheme from converging. The numerical solution 
has the same shape than in the case e > 0: a classical shock followed by a nonclassical shock, 
but this time they have the same speed, thus they remain at the same position. The spike 
is only two cells wide when the classical shocks are reconstructed, and hence not diffused. 
If they are not reconstructed, the spurious classical shock is diffused. 

In conclusion, this test case shows a limitation of scheme, namely its incapability to 
approach exactly classical shock in which w changes sign, but also show its ability to capture 
finely nonclassical solutions, the gap between the shocks being very thin here. In particular, 
the apparition of the intermediate state is immediate, while it is linked to the ratio of the 
width of the gap at time T and of the cell size Ax when using a random sampling based 
method like the Glimm scheme. 

Test 4: Apparition of nonclassical waves from a smooth initial data 

We now focus on the case of a smooth initial data 

f x‘^(x) = 3 sin(27rx), 

[ w^(x) = 1 + 3cos(87rx), 
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Figure 7: Perturbations of an isolated classical shock, with, from top to bottom, e = 0, 
e = 0.05 and e = 0.1. 
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with (5 = 0.95 and m = 1 and periodic boundary condition. Nonclassical shocks appear 
around time t = 0.011 and then propagate in the solution. On Figures [H [9] and [TOl we 
compare the solution given by the reconstruction schemes (with or without reconstruction 
of classical shocks) at time t = 0.015, t = 0.06 and t = 1. The space interval contains 1000 
cells per unit interval and CFL number is set to 0.45. The reference solution is given by the 
Glimm scheme with 8 000 cells. We see that the reconstruction schemes capture accurately 
the nonclassical shocks. The result are poorer in smooth areas, because in those regions 
the reconstruction schemes behave as the Lax-Friedrichs scheme which is quite diffusive. 
This can be improved by using another scheme on interfaces where no reconstruction is 
performed. The use of a moving mesh encourages to chose a central scheme like the Nessyahu 
and Tadmor scheme |NT90| . This scheme is second order accurate in smooth regions of the 
solutions. It is both easy to implement and fast, as it does not use any information on the 
Riemann problems. 



Figure 8: Solution of test 4 at time t = 0.015. 


Test 5: Long time simulation in the maximal dissipative case 

We reproduce here the test case of [CL03] for which /3 = 1 and the final time is large. The 
case j3 = 1 corresponds to the case where the entropy dissipation is zero through nonclassical 
shocks. It is a limit case in the theoretical framework of |LeF02| . As the Riemann solver is 
perfectly defined for /3 = 1, it can be explored numerically. 
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Figure 9: Solution of test 4 at time t = 0.06. 



Figure 10: Solution of test 4 at time t = 0.1. 
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The initial data is 1-periodic with 




'(0.3,0.4) 

< (0.15,-0.2) 
.(0.1,0.4) 


if 0 < X < 0.3, 

if 0.3 < X < 0.3+ 2/3, 

if 0.3+ 2/3 < X < 1. 


Its mean value is null. The parameter m is fixed to 0.05. On Figure [TT] we plot the solution 
at time t = 40 with 2 000 and 8 000 points per unit interval. We can see that at that time 
w changes sign three times, instead of two in the initial solution, and does not converge 
to zero. The position of the nonclassical shocks are different with the two schemes. Let 





Figure 11: Solution of test 5 at time t = 40 (the shape of v is similar) 

us recall that the Glimm scheme is based on a random sampling of the solution, thus two 
realizations of the same test give different results. On Figure [T21 we plot the histogram of 
the hrst nonclassical shock position at time t = 20 for 100 independent realizations of the 
Glimm scheme with 2 000 cells (bottom) and the comparison of the two schemes at the same 
time (with 8 000 points). Moreover, the structure with 2 nonclassical shocks very close to 
each other only appears in 31 out of those 100 realizations. Its indicates that the width 
of this structure is small (and probably smaller than the one appearing in the realization 
of [CL03] L 
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Figure 12: Solution of test 5 at time t = 40 (the shape of v is similar) 
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